## Libraries
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ ggplot2 3.3.0     ✓ purrr   0.3.4
## ✓ tibble  3.0.1     ✓ dplyr   0.8.5
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
# library(purrr)
# library(httr)
# library(magrittr)

library(raster)
## Loading required package: sp
## 
## Attaching package: 'raster'
## The following object is masked from 'package:dplyr':
## 
##     select
## The following object is masked from 'package:tidyr':
## 
##     extract
library(rgdal)
## rgdal: version: 1.4-8, (SVN revision 845)
##  Geospatial Data Abstraction Library extensions to R successfully loaded
##  Loaded GDAL runtime: GDAL 3.0.4, released 2020/01/28
##  Path to GDAL shared files: 
##  GDAL binary built with GEOS: TRUE 
##  Loaded PROJ.4 runtime: Rel. 7.0.0, March 1st, 2020, [PJ_VERSION: 700]
##  Path to PROJ.4 shared files: (autodetected)
##  Linking to sp version: 1.4-1
library(lidR)

library(functional)
library(RColorBrewer)

library(glue)
## 
## Attaching package: 'glue'
## The following object is masked from 'package:raster':
## 
##     trim
## The following object is masked from 'package:dplyr':
## 
##     collapse
library(furrr)
## Loading required package: future
## 
## Attaching package: 'future'
## The following object is masked from 'package:raster':
## 
##     values
# future::plan(multiprocess)
library(sf)
## Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 7.0.0
crs <- st_crs(32650)
crs$proj4string
## [1] "+proj=utm +zone=50 +datum=WGS84 +units=m +no_defs"
years <- c(2013,2014)
ws <-c(18)
dir_data <-"../Data/lidar"

dir_chm <- paste(dir_data,"raster",sep="/")
dir_treetops <- paste(dir_data,"vector","treetops",years,sep="/")

dir_crowns <- paste(dir_data,"raster","crowns",years,sep = "/")
walk(dir_crowns,dir.create)
## Warning in .f(.x[[i]], ...): '../Data/lidar/raster/crowns/2013' already exists
## Warning in .f(.x[[i]], ...): '../Data/lidar/raster/crowns/2014' already exists
filename_crowns_index <- paste(dir_data,"raster","crowns","index.txt",sep="/")
filename_crowns_errors <- paste(dir_data,"raster","crowns","errors.txt",sep="/")
## Read rasters
chm_filenames <- paste0(dir_chm,"/raster",years,".tif")
chms <- map(chm_filenames,~raster(.,crs=crs$proj4string))
read_treetops <- function(year,ws){
  filename <-paste0("treetops_lmf_ws",ws,".json")
  filename_full <- paste(dir_data,"vector","treetops",year,filename,sep = "/")
  readOGR(filename_full,p4s = crs$proj4string)
}
treetops <-
  map(years,~read_treetops(.,ws))
## OGR data source with driver: GeoJSON 
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/treetops/2013/treetops_lmf_ws18.json", layer: "2013_ws18"
## with 2510 features
## It has 2 fields
## Warning in readOGR(filename_full, p4s = crs$proj4string): p4s= argument given as: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
##  and read as: +proj=longlat +datum=WGS84 +no_defs
##  read string overridden by given p4s= argument value
## OGR data source with driver: GeoJSON 
## Source: "/home/lune/Data/ai4er/mres/lidar/vector/treetops/2014/treetops_lmf_ws18.json", layer: "2014_ws18"
## with 2583 features
## It has 2 fields
## Warning in readOGR(filename_full, p4s = crs$proj4string): p4s= argument given as: +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs
##  and read as: +proj=longlat +datum=WGS84 +no_defs
##  read string overridden by given p4s= argument value
# # subset tall trees
treetops <-
  treetops %>%
  # map(~subset(.,.$Z > 60))%>%
  map(st_as_sf)%>%
  {.}
treetops
## [[1]]
## Simple feature collection with 2510 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 587149.8 ymin: 547214.2 xmax: 588598.8 ymax: 548713.2
## CRS:            +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
##    treeID        Z                  geometry
## 1       1 50.83929 POINT (587182.8 548713.2)
## 2       2 66.56500 POINT (587307.8 548713.2)
## 3       3 65.89000 POINT (587513.2 548713.2)
## 4       4 46.18609 POINT (587549.2 548713.2)
## 5       5 66.71000 POINT (587596.2 548713.2)
## 6       6 62.41000 POINT (587615.8 548713.2)
## 7       7 45.32000 POINT (587692.8 548713.2)
## 8       8 60.35218 POINT (587780.8 548713.2)
## 9       9 51.00000 POINT (587933.2 548713.2)
## 10     10 60.14000 POINT (588070.2 548713.2)
## 
## [[2]]
## Simple feature collection with 2583 features and 2 fields
## geometry type:  POINT
## dimension:      XY
## bbox:           xmin: 587149.2 ymin: 547213.2 xmax: 588599.2 ymax: 548713.8
## CRS:            +proj=utm +zone=50 +datum=WGS84 +units=m +no_defs +ellps=WGS84 +towgs84=0,0,0
## First 10 features:
##    treeID        Z                  geometry
## 1       1 50.31984 POINT (587182.8 548713.8)
## 2       2 51.99703 POINT (587211.2 548713.8)
## 3       3 61.81390 POINT (587394.2 548713.8)
## 4       4 64.12300 POINT (587478.8 548713.8)
## 5       5 64.88100 POINT (587514.2 548713.8)
## 6       6 63.59278 POINT (587601.8 548713.8)
## 7       7 37.66848 POINT (587678.8 548713.8)
## 8       8 37.49701 POINT (587692.8 548713.8)
## 9       9 59.36400 POINT (587753.8 548713.8)
## 10     10 50.78204 POINT (587933.8 548713.8)
# plot(chms[[1]])
# plot(treetops[[1]]$geometry,add=T)
# plot(chms[[2]])
# plot(treetops[[2]]$geometry,add=T)
plot(chms[[1]])
plot(treetops[[1]]$geometry,add=T)
plot(treetops[[2]]$geometry,col="red",pch="+",add=T)

dm <- st_distance(treetops[[1]],treetops[[2]])
dim(dm)
## [1] 2510 2583
a <-
  map(array_branch(dm,1), sort) %>%
  as_vector() %>%
  array(dim=rev(dim(dm))) %>%
  t() %>%
  # order()%>%
  {.}
b <-a[order(a[,1]),]

plot(b[,1],pch=".")
points(b[,2],pch='.',col='red')
title(paste("Distance to two closest trees, ordered by 1st, size",ws))

b <-a[order(a[,2]),]

plot(b[,1],pch=".")
points(b[,2],pch='.',col='red')

title(paste("Distance to two closest trees, ordered by 2nd, size",ws))